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Abstract 

We present a detailed description of our calculation of next-to-leading order QCD 
corrections to heavy quark production in e + e~ collisions including mass effects. 
In particular, we study the observables and in the E, EM, JADE and 
DURHAM jet-clustering algorithms and show how one can use these observables to 
obtain m^mz) from data at the Z peak. 
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1 Introduction 



The origin of fermion masses is one of the unresolved puzzles in present high 
energy physics. To be able to answer this question one needs to know the 
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values of the quark masses with high accuracy. However, due to confinement, 
quarks do not appear as free particles in nature, and therefore, the definition of 
their mass is ambiguous. Quark masses can be understood more easily if they 
are interpreted like coupling constants rather than fixed inertial parameters 
and, therefore, quark masses can run if measured at different scales. Moreover, 
in the standard model (SM) all fermion masses come from Yukawa couplings 
and those also run with the energy. This has very important phenomenological 
consequences in Higgs boson searches since the parameter that governs the 
decay rate of the Higgs boson to bottom quarks [1-3] is the running mass of 
the 6-quark at the m# scale, m&(m#), which is much smaller than the quark 
mass extracted at threshold. On the other hand, from the model building 
point of view, to test fermion mass models one has to run masses extracted at 
quite different scales up to the same scale and compare them with the same 
"ruler" [4]. This way, for instance, one can check that, in some grand unified 
models, although the 6-quark mass and the tau mass are very different at 
threshold energies, they could be equal at the unification scale [5,6]. 

However, although predicted by quantum field theory, the possibility of testing 
experimentally the running of fermion masses has not been considered until 
very recently [6-14]. The reason being that for energies much higher than 
the fermion mass the mass effects become negligible for most observables. In 
principle, it is clear that the running of quark masses will be checked, once 
the Higgs is discovered, by measuring directly the Yukawa couplings in Higgs 
decays. We think it is also interesting to explore the possibility of measuring 
quark masses far away from threshold in order to check the running. 

In [7,15-18] it was shown that mass effects in three-jet production at the Z 
peak are large enough to be measured but in [7] we also showed that a next-to- 
leading order (NLO) calculation of three-jet ratios including mass effects was 
necessary in order to extract a meaningful value for the 6-quark mass. This 
is because the leading order (LO) calculation does not distinguish among the 
different quark mass definitions (e.g. perturbative pole mass or MS running 
mass). 

Although jet production of heavy quarks has been considered in a large vari- 
ety of processes, there are very few NLO calculations of jet production taking 
into account complete mass effects. In [19] those were calculated for gluon- 
gluon fusion in proton-antiproton collisions, while in [20] the complete NLO 
corrections were computed for virtual-photon production of heavy quarks in 
deep-inelastic electron-proton scattering. However, until very recently, no full 
calculation of jet production of heavy quarks in e + e~ collisions was available 
at NLO0. In [8,12] we have presented the final results of such a calculation 



Completely inclusive quantities as the total cross-section were fully known at 
order a s [21-24] and some leading quark mass effects in those quantities were also 
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and have shown how can one use the DURHAM (kx) clustering algorithm to 
extract the 6-quark mass at the mz scale. Finally in [14] the DELPHI collab- 
oration has presented the measurement of m h {m z )- The obtained value is in 
good agreement with low energy measurements and running of the mass from 
/i — mi, to /i — mz as predicted by QCD with 5 flavours. Our NLO calculation 
was also used in [14] to check the universality of the strong interactions. In 
the meanwhile two more groups, ref. [36,37] and ref. [38-41], have also pre- 
sented NLO results of jet production of heavy quarks in e + e~ collisions by 
using different calculational schemes. On the other hand also the different ex- 
perimental collaborations have started to use NLO calculations of mass effects 
in order to check the universality of the strong interactions [42-44] and have 
also shown that mass effects are absolutely necessary to explain the data. 

In this paper we present a detailed description of our calculation of the NLO 
corrections to the rates of jet production of massive quarks in e + e~ collisions 
for the different jet-clustering algorithms discussed in [7]. In [14,12] only the 
DURHAM algorithm was used because it has a very good behaviour from 
both the experimental and the theoretical points of view. It is important to 
stress that not all the observables are equally good to study mass effects. In 
particular observables which suffer from huge NLO corrections cannot be used 
to extract any reasonable value of the 6-quark mass. 

Although, in principle, we will concentrate on 6-quark production at the Z- 
peak, some of the results can also be applied to t-t production at the NLC 
(Next-Linear-Collider) . 

In section 2 we review the leading order results, give LO amplitudes in D- 
dimensions and introduce the jet-clustering algorithms we use. In section 3 the 
strategy of the NLO calculation in the framework of the phase-space slicing 
method [45-51] is explained. Section 4 gives all the contributions coming from 
three partons, while in section 5 we give all four-parton contributions and show 
the cancellation of the infrared (IR) divergences among three-parton and four- 
parton contributions. Finally in section 6 we present numerical results for the 
finite parts and review applications of this calculation. In particular we study 
two observables, Rf (the ratio of the three-jet rate containing the 6-quark to 
the three-jet rate of light quarks) and (the ratio of the differential two-jet 
rate containing the 6-quark to the differential rate of light quarks), computed 
in the different jet-clustering algorithms and show how they can be used to 
extract m^mz) from data at the Z-peak. In the appendices we collect some 
of the formulae used in the main calculation. Thus, in appendix A we give the 
most relevant one-loop scalar functions that appear in the calculation of the 
virtual corrections. In appendix B we review the calculation of phase-space 



known up to order c? a [25,26]. Quark mass effects in three-jet final states, however, 
were only known at leading order [24,27-35]. 
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in D-dimensions. Finally in appendix C we collect some of the IR divergent 
integrals which appear in the calculation of four-parton contributions. 



2 The decay Z — > 3jets with heavy quarks at the leading order 

2.1 Decay into two and three jets and jet- clustering algorithms 

The lowest order contribution to the Z-boson decay into three jets in the 
parton picture is given by the decay Z — > bbg. This partonic process has 
an infrared singularity because a massless gluon could be radiated with zero 
energy^. It is well known [52-54], however, that this divergence is exactly 
canceled by the IR-divergent part of the two-parton decay width of the Z- 
boson, Z — > bb, at the order 0(a s ). In the last process the IR divergence 
occurs due to massless gluons running in the loops. 

Thus, in order to define IR-fmite observables, we need to introduce some 
resolution parameter, y c , and to split the full three-parton phase-space into 
the two-jet region, containing soft gluon emission, and the remaining three- 
jet region. The sum of the three-parton decay probability integrated over the 
two-jet part of the phase-space and the two-parton decay width at the order 
0(a s ) defines the IR finite two-jet decay width, Y b 2 -, at the next-to-leading 
order. The integration over the rest of the three-parton phase-space gives the 
three- jet decay width, T\-, at the leading order and it is also IR finite. 

Although both IR-finite T\j(y c ) and r^(?/ c ) depend on the resolution param- 
eter, it is obvious that the sum of the two-jet and three-jet decay widths is 
independent of y c and is given by the total (inclusive) decay width of the 
Z-boson into the heavy quarks at the next-to-leading order: 

r b = r b 2j (y c ) + I* (y c ) . 

Therefore, at order 0(a s ), knowing T b and r^(?/ c ) we can obtain Y b 2; j{y c ) as 
well. The analytical expression for the T b with quark-mass corrections can be 
found in [7,25,26,55-57]. 

In the same way, at the next-to-next-to-leading order we can write for the 



4 In the limit of massless quarks there are collinear singularities as well. In the 
massive case the quark-gluon collinear singularities are softened into logarithms of 
the quark mass. Gluon-gluon collinear singularities appear in our calculations at 
the NLO and will be discussed in the following sections. 
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total width: 

T b = T h 2j {y c ) + T%{y c )+T%{y c ) . (2.1) 

Now the four-jet decay width, T\Ay c ), receives contributions from four-parton 
processes. The three-jet decay width, T\Jy c ), gets contributions from order 
0(a 2 s ) one-loop corrected three-parton processes and the soft part of the four- 
parton processes such that the sum is IR-fmite. Thus, knowing T b , Y\j{y^) and 
F\j(y c ) we can obtain Y b 2 ^{y^) at order 0(a 2 s ). 

Let us discuss now the definition of jets in more detail. The most popular jet 
definitions used in the analysis of the e + e~-annihilation during the last years 
are based on the so-called jet- clustering algorithms. These algorithms have 
to be applied to define jets in both, the theoretical calculations at the parton 
level, and in the analysis of the bunch of real particles observed at experiment. 

In the jet-clustering algorithms jets are usually defined as follows: starting 
from a bunch of n particles^] with four- momenta Pi, (i = l..n) one computes 
a resolution function depending on the momenta of two particles, for example, 

E-E ■ 

y i:j = 2-^(1 - cos%) , 

for all pairs (i, j) of particles. Then one takes the minimum of all and, if it 
satisfies that it is smaller than a given quantity y c (the resolution parameter, 
y-cut), the two particles which define the minimal yij are regarded as belonging 
to the same jetQ. Then they are recombined into a new pseudo-particle with 
the four-momentum defined according to some rule, for example, 

Pk=Pi+ Pj ■ 

After this first step one has a bunch of n — 1 (pseudo)particles and the al- 
gorithm is applied again and again until all the remaining (pseudo)particles 
satisfy the criteria y^ > y c . The final number of (pseudo)particles is the num- 
ber of jets in the considered event and the momenta of the (pseudo)particles 
define the event kinematics. 

In theoretical calculations one can define a jet cross-section or a decay width 
into jets as a function of y c , which are computed at the parton level, by fol- 
lowing exactly the algorithm described above. Jet-clustering algorithms lead 
automatically to IR finite quantities. Furthermore, it has been shown in the 

5 In what follows we will use the word "particles" for both partons and real particles. 

6 The assignment of particles to the jet in the recently proposed Cambridge scheme 
is more involved, see refs [58-60] for details. 
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Table 1 



The jet-clustering algorithms 



Algorithm 


Resolution function, 


Combination rule 


EM 


2 (Pi -Pj)/ S 


Pk = Pi + Pj 


JADE 


2E i E j {\ - cos%)/s 


Pk=Pi+ Pj 


E 


{Pi+Pj?/s 


Pk=Pi+ Pj 


DURHAM 


2mm{E?,E])(l -cos%)/s 


Pk=Pi+ Pj 



literature, using Monte Carlo models, that, for some of the algorithms, the 
passage from partons to hadrons (hadronization) does not change much the 
behaviour of the jet-observables [61], thus allowing to compare theoretical pre- 
dictions with experimental results. We refer to [61,62] for a detailed discussion 
and comparison of different successful jet-clustering algorithms used in e + e~ 
annihilation in the case of massless quarks. Hadronization corrections for the 
case of massive quarks have also been studied for the DURHAM algorithm [14]. 

In this paper we will present results and compare them for the four jet- 
clustering algorithms listed in table 1, where ^/s is the total center of mass 
energy. Already in ref. [7] we showed that the ^-algorithm has a very peculiar 
behaviour when massive partons are involved. It is because for the same values 
of particle momenta, the resolution function, y^, is significantly shifted, com- 
paring with the massless case, when the quark masses are included. To avoid 
these problems in the massive case we introduced in ref. [7] the iJM-algorithm, 
which is similar to the ^-algorithm. As mentioned above, all observables have 
a completely different dependence on y c for the ^-algorithm and this can serve 
as a good test of both calculations and data analyses. The application of the 
new interesting Cambridge algorithm [58-60] in the massive case is presently 
under study [63-67]. 



It is convenient to parameterize the final IR-fmite result for the three-jet decay 
width of the Z-boson in the following general form (D = 4): 



T b _ 



a 



a. 



gy(b)H v (y c , r b ) + g 2 A (b)H A (y c , r b ) 



(2.2) 



where the functions H V (A)(yc, r b) depend, obviously, on what jet-clustering 
algorithm has been used to define the three-jet region of phase-space. We 
have introduced the notation r b = M 2 /s, where s = m| for the Z-decay. 
Here and in the following M b will stand for the perturbative pole mass of the 
6-quark, while the running mass will be denoted as m b (fi). The vector and 
axial-vector neutral-current couplings gv(f) and of a fermion, /, in the 

Standard Model are given by 
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gv(f) = 2T/ - AQ f s 2 w , g A (f) = -22?, (2.3) 

being sw (cw) the sine (cosine) of the weak mixing angle, T| the third com- 
ponent of the weak isospin and Qf the electric charge of the fermion. 



We can expand the H V (a) functions in the strong coupling constant, a s 

H V (A) (y c , n) = H { y ] (A) (y c , r b ) + —H$ A) (y c , r h ) H , (2.4) 



where Hy^ and Hy) A s are the LO and NLO contributions, respectively. Fur- 
thermore, in order to see more clearly the size of quark-mass effects, when 
the masses are small with respect to the center of mass energy, we found it 
convenient to rewrite these functions factoring the dominant r^-dependence 

H® A) (y c ,r b ) = A v \ A) (y c )+r b B$ A) (y c ,r b ) . (2.5) 



Note that eq. (2.5) is not an expansion in r b because the functions Byl A , (y c , r b ) 
contain the exact residual ^-dependence. 

In the limit of zero quark masses, r b — 0, chirality is conserved and the two 
lowest-order functions Hy\y c ,r b ) and H A \y c ,r b ) become identical 

H v 1 \y c ,0) = H A °\y c ,0) = AW(y c ) . (2.6) 



This is not true at the NLO where the vector and the axial-vector functions, 
Ay and A A \ are not equal anymore due to the small contribution of the 
triangle diagram [68] V12 in fig. 2. 

The LO and NLO massless contributions, the A^ and Ay functions^, were 
calculated in [61] using the results from [69]. In [7] the lowest order mass 
effects, functions By^ } have been calculated analytically for the EM jet- 
clustering algorithm and numerically for JADE, DURHAM and E algorithms. 
For completeness, we review in the next section the LO results. Results for the 
NLO heavy quark contribution, in particular the functions By) A -, in eq. (2.5) 
will be presented later. 

Concluding this section we would like to make the following remark. In this 
paper we discuss the Z-boson decay. In LEP experiments one studies the 
process e + e~ — > {Z r f) — > bb + ■ ■ ■ and, apart from the resonant Z-exchange 
cross-section, there are contributions from the pure 7-exchange and from the 

7 Note that with our choice of the normalization A^\y c ) = ^A(y c ) and Ay\y c ) = 
\B(y c ) where A(y c ) and B(y c ) are defined in [61]. 



7 



7 — ^-interference. The non-resonant 7-exchange contribution at the peak is 
less than 1% for muon production and in the case of 6-quark production there 
is an additional suppression factor Q\ = 1/9. In the vicinity of the Z-peak 
the interference is also suppressed because it is proportional to Q b (s — m|). 
However, away from the Z resonance these contributions are important and 
should be properly taken into account. 

The extension of our calculation for the cross-section of e + e~ annihilation into 
three jets can be done [55] as follows 

^3j(s,y c ,r b ) = — \a v (s)H v (y c ,r b ) + (T b A (s)H A (y c ,r b ) 

J n L 



with the same Hy^ijjc^'b) functions introduced in eq. (2.2) for the Z-decay. 
The functions Cy(A)( s ) have the form (D — 4): 



4(*) 

where 



47ra 2 [ n2 ^ 2 1 



+ 



QtQt + 2Q e Q b 2 2 g v (e)g v (b)Re(x(s)) 

1 



<7v(&) 2 ^(e) 2 + ^(e) 2 



47Ttt 2 



s \ 64s^c^ 



^(6) 2 ^y(e) 2 + ^(e) 2 



(2.7) 



X(s) 



s — m| + im z Y z 



Also, the numerically important QED initial-state radiation (ISR) should be 
taken into account in the realistic analysis of the experimental data. The 
cross-section of the three-jet production including ISR can be written as a 
convolution 



v b 3j(s,y c ,r b ) = J a b 3j (s' ,y c ,r b )F ISR (s' / s)ds , 



(2.8) 



where F ISR (s'/s) is the well-known QED radiator [70] for the total cross- 
section and s' denotes the invariant mass of the jet final state. 



2.2 Z — > 3 jets at tree level 



In the NLO calculation of the Z-decay width into three-jets, which will be 
presented in the next sections, we encounter ultraviolet (UV) and infrared 
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singularities. We use dimensional regularization to regularize both types of 
divergences and therefore most of the calculation has to be done in arbitrary 
D-dimensions. The calculation of Y b ^{y c ) at order a s is a pure tree- level cal- 
culation and it does not have any IR problem, since the soft gluon region is 
excluded from the three-parton phase-space. Therefore, the calculation can be 
safely done in four dimensions. However, the later use of the LO results in some 
steps of the NLO calculation require the tree-level results in D-dimensions that 
we present in the following. 




Ta Tb 



Fig. 1. Tree-level diagrams contributing to Z — > bbg. 

In D = 4 — 2e dimensions the Born transition probability for the three-parton 
process, Z — > b b g, summed over final colours and spins and averaged over 
initial spins is equal to 

|M * I = 16^3 = 2^' (2 ' 9) 



where g the SU{2) coupling constant, g s is the strong coupling constant, Cf = 
4/3 and N c = 3 are SU(3) group invariant factors and 



Ttbg(yia,y2a) = 16 



2/132/23 



T bi + gUb)(l-e) 



i (1 _ c) (m + y*)- 



V2/23 Vl3, 



+9 2 A (b)(l-e) 



2 V 2/23 2/13 / 



(2.10) 



with 



2/13 = 



2(pi-p 3 ) 



2/23 



2(^2 -Ps) 



where p±, p 2 and p% are the four-momenta of the quark, the antiquark and the 
gluon, respectively. The function ^(^13,3/23) defines the three-body phase- 
space boundaries (See eq. (B.3)) and the dimensionless function 

T bi = <&(&)(! - e + 2r b ) + g 2 A (b)(l - e)(l - 4r b ) , (2.11) 
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is related to the Born two-parton process, Z — > b b, transition probability 
through 



9 2 N ( 



(2.12) 



W w 3 - 2e 



Averaging over the Z-boson polarizations we took into account that in arbi- 
trary space-time dimensions the number of spin degrees of freedom is D — 1 = 
3-2e. 

In the limit of zero energy of the emitted gluon, both invariants yi% and 2/23 
approach zero and the first term in the eq. (2.10) develops an IR singularity. 
The IR-finite three-jet decay width of the Z-boson at LO is given by the 
integral of the transition probability T bbg (yi 3 , 1/23), eq. (2.10), over the three- 
jet region of the phase-space, where the soft-gluon region is excluded: 



r fe 



1 



JdPS(3)J2\M bh \ 2 e c (y c ;y 13 ,y 23 ) 



(2.13) 



=0, s=mi. 



The function 9 c (y c ] y 13 , 1/23) introduces the appropriate cuts for each of the 
jet-clustering algorithms and defines the three-jet region. 

At the lowest order, E, JADE and EM give the same three-jet rates for mass- 
less particles, because in this case there is no parton recombination involved 
and all the three schemes have the same resolution parameter. However, al- 
ready at order a 2 they give different results since after the first recombination 
the pseudo-particles are not massless anymore and the resolution functions 
are different. For massive quarks the three algorithms, E, JADE and EM are 
already different at order a s . The DURHAM algorithm is, of course, com- 
pletely different from the other algorithms we use, both in the massive and 
the massless cases. 

Tables 2, 3, 4 and 5 present the dependence of the function A^(y c ) and the 
ratio of functions By^(y c ,r b )/A^(y c ) on the jet-resolution parameter y c for 
the different jet schemes and for several relevant values of M b . It is clearly 
seen from the tables that both leading order functions B$ and B^ have a 
soft dependence on M b for the region of y-cut, y c > 0.01, which is relevant for 
experimental studies. 
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Table 2 

The leading order functions (see eqs. (2.2, 2.4, 2.5)) for the E algorithm. 



Vc 



0.01 
0.02 
0.03 
0.04 
0.05 
0.06 
0.07 
0.08 
0.09 
0.10 



Vc 



0.01 
0.02 
0.03 
0.04 
0.05 
0.06 
0.07 
0.08 
0.09 
0.10 



A (o) 



18.627 
12.209 
9.069 
7.118 
5.765 
4.764 
3.990 
3.374 
2.873 
2.458 



18.627 
12.209 
9.069 
7.118 
5.765 
4.764 
3.990 
3.374 
2.873 
2.458 



4 0) /A(°) 



M b = 3GeV 



37.811 
22.195 
16.917 
13.540 
11.390 
9.932 
8.889 
8.118 
7.540 
7.108 



M b = 4GeV 



39.353 
21.645 
16.314 
13.560 
11.466 
9.981 
8.920 
8.138 
7.551 
7.113 



M b = 5GeV 



42.432 
21.723 
15.840 
13.135 
11.452 
10.043 
8.960 
8.162 
7.565 
7.118 



M h = 3GeV 



32.506 
17.188 
12.084 
8.838 
6.797 
5.436 
4.478 
3.785 
3.280 
2.914 



M h = 4GeV 



33.827 
16.532 
11.414 
8.812 
6.838 
5.455 
4.486 
3.785 
3.274 
2.905 



M h = 5GeV 



36.578 
16.462 
10.847 
8.321 
6.777 
5.481 
4.495 
3.785 
3.267 
2.893 



3 Three-jet ratios at next-to-leading order: overview of the calcu- 
lation 



Like for many QCD processes it turns out that the leading order calculation of 
the Z-boson decay into three-jets, described in the previous section, does not 
match the experimental precision. The numerical results for the leading order 
calculation were presented in [33,34,7] (see also section 2.2). It is known [61] 
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Table 3 

The leading order functions (see eqs. (2.2, 2.4, 2.5)) for the EM algorithm. 



Vc 




4 0) /a(°) 






M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


18.628 


-18.708 


-18.778 




-18.444 


0.02 


12.209 


-12.406 


-13.178 




-13.529 


0.03 


9.069 


-10.081 


-10.637 




-11.167 


0.04 


7.118 


-9.464 


-9.480 




-9.813 


0.05 


5.765 


-9.162 


-9.161 




-9.168 


0.06 


4.764 


-9.018 


-9.016 




-9.013 


0.07 


3.990 


-8.973 


-8.970 




-8.966 


0.08 


3.374 


-8.997 


-8.993 




-8.988 


0.09 


2.873 


-9.073 


-9.068 




-9.062 


0.10 


2.458 


-9.192 


-9.186 




-9.178 


Vc 




Bf/A^ 






M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


18.628 


-23.675 


-23.678 




-23.261 


0.02 


12.209 


-17.214 


-17.931 




-18.215 


0.03 


9.069 


-14.764 


-15.269 




-15.736 


0.04 


7.118 


-14.042 


-14.008 




-14.279 


0.05 


5.765 


-13.648 


-13.598 




-13.544 


0.06 


4.764 


-13.419 


-13.370 




-13.308 


0.07 


3.990 


-13.295 


-13.247 




-13.184 


0.08 


3.374 


-13.246 


-13.198 




-13.135 


0.09 


2.873 


-13.255 


-13.205 




-13.142 


0.10 


2.458 


-13.310 


-13.260 




-13.195 



that next-to-leading corrections to the decay of the Z-boson into massless 
quarks are very significant numerically. In addition, at the leading order we 
do not have any idea what value of the quark mass should be taken. The 
various perturbative definitions of the quark mass (pole mass, running mass 
etc.) differ only at the next-to-leading order. However, the numerical differ- 
ence when several mass definitions are used in the leading order prediction 
for three-jet observables is large [7]. It is again an indication that next-to- 
leading corrections are important. Therefore, to use consistently the quark 
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Table 4 

The leading order functions (see eqs. (2.2, 2.4, 2.5)) for the JADE algorithm. 



Vc 




4 0) /a(°) 






M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


18.627 


-34.488 


-31.314 




-28.526 


0.02 


12.209 


-25.045 


-24.050 




-22.805 


0.03 


9.069 


-20.158 


-19.968 




-19.527 


0.04 


7.118 


-17.410 


-17.308 




-17.191 


0.05 


5.765 


-15.675 


-15.596 




-15.501 


0.06 


4.764 


-14.491 


-14.437 




-14.384 


0.07 


3.990 


-13.613 


-13.577 




-13.529 


0.08 


3.374 


-12.924 


-12.902 




-12.870 


0.09 


2.873 


-12.393 


-12.373 




-12.346 


0.10 


2.458 


-11.979 


-11.964 




-11.943 


Vc 




Bf/A^ 






M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


18.627 


-39.412 


-36.144 




-33.244 


0.02 


12.209 


-29.823 


-28.753 




-27.415 


0.03 


9.069 


-24.817 


-24.560 




-24.035 


0.04 


7.118 


-21.968 


-21.802 




-21.606 


0.05 


5.765 


-20.142 


-20.002 




-19.832 


0.06 


4.764 


-18.876 


-18.764 




-18.638 


0.07 


3.990 


-17.921 


-17.829 




-17.711 


0.08 


3.374 


-17.161 


-17.086 




-16.986 


0.09 


2.873 


-16.564 


-16.492 




-16.399 


0.10 


2.458 


-16.088 


-16.022 




-15.937 



mass definition, and to provide accurate predictions for three-jet observables, 
the next-to-leading order, 0(a^), calculation for massive quarks has to be 
included. 

The main difficulty of the NLO calculation is the appearance at the interme- 
diate stages, in addition to ultraviolet divergences, of infrared and collinear 
singularities due to massless gluons. The Bloch-Nordsieck and Kinoshita-Lee- 
Nauenberg theorems [52-54] assure, however, that observable jet cross-sections 
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Table 5 

The leading order functions (see eqs. (2.2, 2.4, 2.5)) for the DURHAM algorithm. 



Vc 




bP/a^ 






M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


7.836 


-27.862 


-27.092 




-26.202 


0.02 


5.127 


-20.038 


-19.800 




-19.509 


0.03 


3.815 


-16.760 


-16.652 




-16.517 


0.04 


3.003 


-14.889 


-14.833 




-14.762 


0.05 


2.440 


-13.659 


-13.628 




-13.589 


0.06 


2.023 


-12.780 


-12.764 




-12.743 


0.07 


1.701 


-12.117 


-12.110 




-12.101 


0.08 


1.443 


-11.598 


-11.597 




-11.595 


0.09 


1.233 


-11.180 


-11.183 




-11.186 


0.10 


1.058 


-10.836 


-10.842 




-10.849 


Vc 










M b = 3GeV 


M b = 4GeV 


M b 


= 5GeV 


0.01 


7.836 


-32.216 


-31.337 




-30.316 


0.02 


5.127 


-24.256 


-23.934 




-23.537 


0.03 


3.815 


-20.880 


-20.698 




-20.471 


0.04 


3.003 


-18.930 


-18.807 




-18.652 


0.05 


2.440 


-17.634 


-17.541 




-17.423 


0.06 


2.023 


-16.698 


-16.623 




-16.527 


0.07 


1.701 


-15.985 


-15.922 




-15.841 


0.08 


1.443 


-15.421 


-15.366 




-15.295 


0.09 


1.233 


-14.964 


-14.915 




-14.851 


0.10 


1.058 


-14.584 


-14.539 




-14.481 



are infrared finite and free from collinear divergences. 



At the NLO the three-jet cross-section has contributions from three-parton, 
figs. 2 and 3, and four-parton, figs. 5, 6 and 7, final states. 

The IR and collinear singularities of the NLO one-loop Feynman diagrams 
cancel against divergences that appear when the differential cross-section for 
four-parton production is integrated over the region of phase-space where ei- 
ther one gluon is soft or two gluons are collinear. However, the singularities 
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that appear in the intermediate steps of the calculation should be treated 
properly. We used dimensional regularization to regularize both UV and IR 
divergences [71-73] because this regularization preserves the QCD Ward iden- 
tities. 

The three-parton transition amplitude can be expressed in terms of a few 
scalar one-loop integrals. The result contains poles in e = (4 — D)/2, where 
D is a space-time dimension. Some of the poles have UV origin and other 
correspond to the IR singularities. All UV divergences are removed by appro- 
priate renormalization of the parameters of the QCD Lagrangian (coupling 
constant, mass and wave functions). After UV renormalization and taking the 
interference with the tree-level three-parton amplitude we obtain analytical 
expressions which have a part containing the IR poles (single and double) and 
a finite contribution. The singular part is proportional to the leading order 
transition probability. 

The four-parton transition probabilities are split in two parts. The first one, 
the so-called soft part, contains singularities when one of the massless partons 
in the final state is soft, or two massless final partons are collinear. The second 
part, denoted as hard, is free from any potential singularities. Then, to cancel 
the IR/collinear divergences against those appearing in the three-parton con- 
tribution we work in the context of the phase-space slicing method [45-51]. 
In this method the analytical integration over a thin slice at the border of 
the phase-space, which contains the IR/collinear singularities, is performed in 
/^-dimensions by using approximate expressions for the transition probability 
expanded in the slice-parameter that defines the size of the slice. The result 
is added to the virtual corrections. The sum becomes free of singularities and 
can be integrated numerically for D = 4 in the three-jet region of the three- 
body phase-space (defined by the jet-clustering algorithms considered in the 
previous section: EM, JADE, E and DURHAM), but it depends on the small 
slice-parameter. The integration of the hard part over the three-jet region 
of the four-body phase-space (again defined by the jet-clustering algorithms 
considered above) is done numerically for D = 4. 

The sum of the three-parton and four-parton contributions to the three-jet 
decay width of the Z boson should be independent of the parameter that 
defines the slice as long as it is small enough not to spoil the approximations 
used. The independence of this parameter has been checked in our calculations. 
Finally, we obtain the finite functions H v {i) ci r b ) and HA(y c , r b) in eq. (2.5) at 
order a s . 

Our classification of the transition probabilities is similar to the one used by 
Ellis, Ross and Terrano [69] which is based on a colour factor classification. The 
way the cancellation of IR divergences occurs can be seen by representing the 
different transition probabilities as the different cuts one can perform in the 
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three-loop bubble diagrams contributing to the Z-boson selfenergy. Therefore, 
we assign both three- and four-parton transition probabilities to the different 
classes defined by the six bubble diagrams depicted in fig 4. The advantage 
of this classification is that the cancellation of IR divergences occurs for each 
class of diagrams separately and, therefore, numerical tests can be performed 
independently for the IR finite result obtained in each class. 



4 Three-parton contributions 

4-1 Classification of diagrams 

The complete set of diagrams describing the one-loop radiative corrections to 
the process 

Z(q)^b( Pl )+b(p 2 )+g(p 3 ) , (4.1) 



is shown in figs. 2 and 3. They contribute to the three-jet decay rate at 0{a 2 s ) 
through their interference with the lowest order bremsstrahlung diagrams Ta 
and Tb in fig. 1. 




VI V2 V3 V4 





V9 V10 Vll V12 



Fig. 2. Radiative corrections to the process Z — > bbg. Diagrams VI to V12 contribute 
at 0{a 2 s ) through their interference with the lowest order bremsstrahlung diagrams 
Ta and Tb. 
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SI 



S2 



S3 



S4 




S5 



S6 



S7 



S8 



Fig. 3. Selfenergy diagrams. Graphs involving ghosts and similar in structure to 55 
and 56 have not been shown. 

We denote by Via and Sja the interference of diagrams Vi and Sj with Ta 
averaged (summed) over initial (final) colours and spins, 



Interference with Tb is denoted in the same way by Vib and Sjb. 

As usual in QCD calculations in the MS scheme at LEP energies, we are 
working in a theory with Np = 5 quark flavours in which only the 6-quark is 
considered massive. Therefore, in all closed quark loops in fig. 3 all five quarks 
are running in the loop. The contribution of the triangle diagram V12 in fig. 2 
is peculiar. This contribution is due entirely to the fact that the top-quark 
and the bottom-quark are not degenerated [68] (lighter quarks are effectively 
degenerated), thus, both top and bottom quarks are kept in the loop. 




Class A Class B Class C Class D Class E Class F 

Fig. 4. Bubble classification relating virtual and real contributions at O(a^) accord- 
ing to their divergent structure. Again we do not show diagrams similar to Class D 
with ghosts in the loop. 

It turned out very convenient to classify the different one-loop three-parton as 
well as four-parton contributions based on the six three-loop bubble diagrams 
presented in fig. 4. Different cuts of the same diagram represent either a tree- 
level four-parton transition probability, or a one-loop three-parton transition 
probability, or a two-loop two-parton transition probability. Contributions, 
corresponding to different bubble diagrams are proportional to different colour 
factors. Then, the cancellation of IR divergences occurs separately between 
contributions obtained from the same class of diagrams. 



Via = Y, Re { Vi * Ta+ }> s J a = Y, Re i S j * Ta+ } ■ 



(4.2) 
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The classification of all one-loop contributions to the decay Z — > bbg is shown 
in table 6. In every column of the table we group contributions with the 
same colour factor which is shown above the column. Here Tr = 1/2. The 
transition probabilities in the third and the fourth rows of the table can be 
obtained from the ones in the first and second rows by the interchange of the 
quark and antiquark momenta. 

Table 6 

The interchange table relating the graphs for Z — > bbg. The third and the fourth 
rows are obtained from the first and the second rows by the interchange of the quark 
and the antiquark momenta (p\ <-> P2). All colour factors are given up to a global 
constant Cf Nc- 



Class A 


Class B 


Class C 


Class D 


Class E 


Class F 


C F 


C F - \N C 


N c 


N c 


Tr 


Tr 


Via V3a Sla S3a 


V5a V9a 


V7a Vila 


S5a 


S7a 


V12a 


Vlb V3b Sib S3b 


V5b V9b 


V7b Vllb 


S5b 


S7b 


V12b 


V2a V4a S2a S4a 


V6a VlOa 


V8a Vila 


S6a 


S8a 


V12a 


V2b V4b S2b S4b 


V6b VlOb 


V8b Vllb 


S6b 


S8b 


V12b 



One-loop contributions from Class A and Class B have at most simple IR poles 
for the case of massive quarks. This is because the corresponding four-parton 
contribution to the three-jet final state can have a singularity only due to the 
radiation of a single soft gluon. 



Class C one-loop contributions involving the three-gluon vertex has, in addi- 
tion, a double infrared pole. In this case the two gluons in the corresponding 
diagram for the process Z — > bbgg involving the three-gluon vertex can be 
collinear at the same time that one of them is soft. 

Contributions from Class D and Class E with a light quark running in the 
inner loop of the bubble have only simple IR poles coming from the gluon 
selfenergy. 

Class E with massive quarks in both loops gives a finite result. 

Class F results in only a small IR-finite contribution to the axial function 
(see eq. (2.4)). 

The calculational procedure of the one-loop contribution is as follows. First, 
we perform all Dirac trace calculations for the matrix element squared. Then, 
as usual, the result is expressed as a sum of Lorentz scalars containing several 
loop-integrals with different number of propagators and with scalar products 
of the loop momentum and external momenta in the numerator. We use the 
standard Passarino-Veltman reduction procedure [74] to reduce all of these 
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vector and tensor loop integrals to simpler scalar n-point functions. Then the 
final result can be written analytically in terms of two IR divergent four- 
propagator (box) integrals, five three-propagator integrals and a number of 
simple two-point and one-point scalar integrals. The analytical results for these 
scalar loop integrals are summarized in appendix A. 

As mentioned above, after dimensional regularization, both UV and IR singu- 
larities appear as poles e. The next step is renormalization and cancellation 
of UV divergences. 



4-2 Renormalization counterterms 



Renormalization is accomplished by using a mixed scheme, on-shell for quark 
masses and MS for the gauge coupling. Then, to subtract the UV divergences 
we make the following replacements in the Born decay widths, 



M h ^M h ll-C, 



g s ^g s {l + — 



An 



^N c - ^T R N F 



(4.3) 
(4.4) 



where A = (A7i) € /(eT(l — e)). Thus, the counterterm that should be added to 
the one-loop decay width to eliminate all the UV divergences has the form, 



2r u Q 



„ 2 



-T R N F - ^rN c 



J dPS(3)T big (y 13 ,y 23 ) , (4.5) 



where we defined, 



2^16(^3-26' ^ S r(l-e) ' 1 J 



and T b i g {yi 3l y2z) is the square of the tree-level single bremsstrahlung ampli- 
tude, see eq. (2.10). 

After renormalization has been performed, the three-parton differential tran- 
sition probability has a finite part and divergent terms proportional to the IR 
poles. 
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4-3 Infrared divergent contributions 



The IR divergent piece of the one-loop contribution to the decay width Z 
bbg reads (See appendix A for notation and details): 



r2 ta «i = 2r a2 JdPS(3)T b - bg (y 13 , 



2/23 ) 



x \~Cf --[Cf-^t) - ^_, Q logc 12 
-N c 




+ l T « 



1 1/11 

e 2 e V 6 Via V23 , _ 

(4.7) 



(N F -l)^- + \ogr b 



Note that in eq. (4.7) we also included a term proportional to T R and logr&, 
which is finite for a non-zero quark mass (r b ^ 0). This term comes from 
the gluon self-energy and would appear as an additional pole if the 6-quark 
had been considered massless. It is canceled by a four-parton contribution in 
which a gluon emitted from the initial 6-quarks splits into two heavy quarks. 
It is important to note that similar contributions appear in the case of gluons 
emitted from light quarks. 



5 Four-parton contributions 

5. 1 Classification of Diagrams 

5.1.1 Emission of two real gluons 
The process 

Z(q) - 6(pi) + 6(p 2 ) + g( P3 ) + g(p A ) (5.1) 

is described by the eight diagrams shown in fig. 5 and, therefore, the transition 
probability contains, in principle, 36 contributions. However, many of them are 
related by interchange of momentum labels and, at the end, only 13 transition 
probabilities need to be calculated. 

As before, we denote by the interference of diagram Bi with Bj averaged 
(summed) over initial (final) colours and spins, 

Bij = Y,Re{Bi * Bj + } . (5.2) 
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Bl B2 B3 B4 




Fig. 5. Feynman diagrams contributing to the process Z — ► bbgg. 

In table 7 we give the momentum label interchanges necessary to generate all 
the transition probabilities and we classify them according to the previously 
defined bubble groups. We see from the table that it is sufficient to consider 
the following combinations of transition probabilities, 

Class A bbgg = -Bll + 2521 + 522 + 532 , 

Class B bhgg = -BA1 + 2542 + 562 + 552 , 
Class C bbgg = 2(571 + 572 + 582) , 

Class D bbgg = - (577 + 587) , (5.3) 

plus the interchanges (1 <-> 2), (3 <-> 4) and (1 2) (3 <-> 4). 
Table 7 



The interchange table relating the graphs for Z — > All colour factors are given 
up to global constant Nc- 



label 


Class A 


Class B 


Class C 


Class D 


permutation 


c F 


C F - \N C 


N c 


N c 




Bll B21 B22 B32 


B41 B42 B62 B52 


B71 B72 B82 


B77 B87 


(1-2) 


B44 B64 B66 B65 


B41 B61 B62 B63 


B84 B86 B76 


B88 B87 


(3-4) 


B44 B54 B55 B65 


B41 B51 B53 B52 


B74 B75 B85 


B77 B87 


(1 2) (3 ^ 4) 


Bll B31 B33 B32 


B41 B43 B53 B63 


B81 B83 B73 


B88 B87 



The sum over the two physical polarizations of the produced gluons is accom- 
plished by summing over the polarizations with, 



E^V^-jT , (5.4) 

pol 
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but including in Class D Feynman diagrams like B7 and B8 with "external" 
ghosts in order to take into account the fact that the gluon current is not 
conserved. 

Only gluons attached to external legs (quarks or gluons) can generate infrared 
divergences in the three-jet region. Thus, B32 and B52 are fully finite, B21, 
B22, B42 and B62 are IR divergent only if gluon labeled as 3 is soft while Bll 
and B41 are IR divergent in the three-jet region when any of the gluons labeled 
as 3 or 4 is soft. As commented before, this kind of transition probabilities, 
classified into Classes A and B, contain only simple infrared poles in e since 
for massive quarks the quark- gluon collinear divergences are softened into 
logarithms of the quark mass. On the other hand, double infrared poles appear 
for diagrams of Class C because the gluon-gluon collinear divergences are 
present due to the three-gluon vertex. The situation is similar for B77 and 
B87. These diagrams individually contain double infrared poles. Nevertheless, 
their sum is such that at the end only simple infrared poles survive. The reason 
being that they belong to the same Class D as the diagrams with a selfenergy 
insertion in an external gluon leg. 

5.1.2 Emission of four quarks 
First we consider the process 

Z(q) -> 6(pi) + b(p 2 ) + q(p 3 ) + q(pd , (5.5) 
where q stands for a light quark, which is assumed massless. 




Fig. 6. Feynman diagrams contributing to the process Z — ► bbqq, where q stands for 
a light quark. 

This process is described by diagrams shown in the fig. 6, where the pair 
of bottom-antibottom quarks is emitted from the primary vertex, and two 
similar Feynman diagrams, where the heavy quark pair is radiated off a light 
qq system attached to the Z-boson (the so-called heavy quark gluon splitting 
diagrams) . 

The transition probabilities Cll, C22 and C12, corresponding to the diagrams 
in fig. 6, are assigned to Class E 

Class E bhi! = Cll + C12 + (1 <-> 2) , (5.6) 
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and can be obtained from the bubble diagram, fig. 4, with a heavy quark in the 
external loop and a light quark in the internal loop. Due to the light quarks, 
which are considered massless, the above transition probabilities have soft and 
quark-quark collinear singularities in the three-jet region. In principle, these 
divergences can contain double polesQ. However, as in the case of the soft and 
the gluon-gluon collinear divergences in 577, 588 and 587, all double poles 
cancel in the sum, because these transition probabilities are related to the 
one-loop three-parton contributions with a external gluon self-energy insertion 
(they correspond to different cuts of the same bubble diagram) , which has only 
a simple IR pole. 

The remaining contributions to the decay Z — > bbqq correspond to the heavy 
quark gluon splitting diagrams and to their interference with the diagrams 
shown in fig. 6. They can be understood following the bubble representation 
discussed in the previous section and depicted in fig. 4. 

The contribution of the heavy quark gluon splitting diagrams is IR finite. It is 
obtained from Class E bubble diagrams in fig. 4 with a light quark (external 
loop) and a heavy quark (internal loop). This four-parton contribution con- 
tains numerically large collinear logarithms of the heavy quark mass. However, 
they all cancel if the four-parton part is summed together with the one-loop 
three-parton contributions with light quarks coupled to the Z boson and a 
gluon self-energy insertion with a heavy quark in the loop. These three-parton 
contributions are obtained from the same Class E bubble diagram by cutting 
the gluon and the light quark lines. In principle, one can assign the heavy 
quark quark gluon splitting contributions to either the light quark three-jet 
decay width or to the heavy quark three-jet decay width. Because the can- 
cellations discussed above we choose to include them in the three-jet Z-decay 
width into light quarks. In this way the limit Mb — > can be taken and re- 
sults can be compared with known massless calculations [62,61,69]. Obviously, 
these heavy quark gluon splitting contributions should be subtracted from the 
experimental data before comparing with our theoretical result. 

Finally, we have to consider the contributions corresponding to Class F bub- 
ble diagrams in fig. 4 which are also IR finite. By cutting one gluon line and 
one quark loop in this bubble diagram one obtains one-loop three-parton con- 
tributions arising from the triangle diagrams V12, which, as commented in 
the previous section, produce a tiny difference, even for massless quarks, be- 
tween HAiUci^b) and Hy(y c ,ri,). Cutting the two quark loops in Class F bub- 

8 For exactly massless quarks. It is also possible to regularize these infrared diver- 
gences with a small quark mass. In this case, infrared divergences are softened into 
mass singularities and lead to large logarithms in the quark mass, log(m g //x). Sim- 
ilarly, infrared gluon divergences can be regulated at lowest order by giving a small 
mass, A, to the gluons. At next-to-leading order we would violate gauge invariance 
at the three-gluon vertex. 
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ble diagrams with one heavy and one light quark loops gives the four-parton 
contribution from the interference of diagrams in which a heavy quark pair is 
radiated off a light quark pair with diagrams in which a light quark is radiated 
off a heavy quark pair. This contribution has the same type of cancellation 
that occur for the corresponding three-parton contribution. Numerically it is 
even smaller than the one produced by triangle diagrams. Therefore we have 
neglected it. 

A different approach to avoid the large logarithms has been followed in [36,37], 
where a double 6-tagging is imposed, sacrificing, however, the experimental 
statistics. 




Fig. 7. Feynman diagrams contributing to the decay width of the Z-boson into four 
massive bottom quarks. 

We consider now the decay into four bottom quarks. 



Z(q) - &(pi) + b(p 2 ) + b(p 3 ) + 6(p 4 ) , (5.7) 

The corresponding diagrams are shown in fig. 7. As in the case of the emission 
of two real gluons, from the eight diagrams shown in fig. 7 we should compute 
only twelve different terms in the transition probability out of the 36 terms 
possible as the other terms are related by interchange of momentum labels. 

Thus, it is sufficient to consider the following combinations: 



-(D15 + D28 + D17 + D26 + 2(D18 + D25)) , 
D11 + 2D12 + D22 , 

D13 + 2£>23 + D24 , (5.8) 

plus the interchanges (1 <-> 3), (2 <-> 4) and (1 <-> 3) (2 <-> 4). Due to Fermi 
statistics there is a relative minus sign for diagrams D5 to D8 that is reflected 
in the transition probabilities of Class B bbbb and ensures they vanish when both 
fermions (antifermions) have identical quantum numbers. 



Class B bbbb = 
Class E bibh = 
Class F bbbb = 
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Table 8 

The interchange table relating the graphs for Z — > bbbb. A global Cf Nq factor has 
been factorized from the colour factor. 



label 
permutation 


Class B 
C F - \N C 


Class E 
Tr 


Class F 
Tr 


(1-3) 
(2-4) 
(1 3) (2 4) 


D18 D25 D15 D28 D17 D26 
D45 D16 D15 D46 D35 D62 
D27 D38 D37 D28 D17 D48 
D36 D47 D37 D46 D35 D48 


Dll D12 D22 
D55 D56 D66 
D77 D78 D88 
D33 D34 D44 


D13 D23 D24 
D57 D67 D68 
D57 D58 D68 
D13 D14 D24 



The transition probabilities of Class F are called in the literature [75] singlet 
contributions because they contain two different fermion loops and hence can 
be split into two parts by cutting gluon lines only. The first contribution to 
the vectorial part arises at 0(al) as a consequence of the non-abelian gen- 
eralization of Furry's theorem. Singlet contributions to the axial part appear 
already at 0(a 2 s ). 



5.2 Infrared divergent contributions 



The part of the four-parton decay width which is singular when a gluon (with 
momentum p 3 ) attached to a massive quark becomes soft is given by, 



= 2r° gt -j- J dPS(3)(y u , y 2 dT bh (yu, 1/24) 



2(2tt)°- 1 J 3 3 \V 2 / 2/13 2/23 2/i 2 3 1 3 ') 

+ (1^2) + (3^4) + (1^2)(3^4) . (5.9) 



The part of the four-parton decay, corresponding to the diagrams with a gluon 
splitting into two gluons or into two light quarks, which is singular when one of 
the two gluons (or light quarks) is soft, or two massless partons are collinear, 
is given by (See appendices B and C for notation and details), 



r5 s ; = 2rV 2 J dPS(3)(y lu , 2/234) 

1 7T 

yu~ e dyu J dv{v{l-v))- e — j dffwT^ff 



w 

x I 



25 



X 



rji / x ) Nc ( 1/134 . 1/234 . 1/134 . 1/234 , ,, .. . 

T b i g {yizA, 1/234) <j — ( — + — + — + - 4 + 2u(l - u) 



I 2! V yia 

+T i? (iV F - 1) 



1/23 



1/14 V2A 



V 2 + (l-v) 



+32 J^— Tbl v{\ - v) (2 cos 2 0' - -1—) 



x(^(i-c)-r fl (Jv F -i)f + o(i^) 



(5.10) 



After integration over the angle the last term, proportional to the lowest or- 
der transition probability T bb , vanishes. The other terms become, as expected, 
proportional to the well-known Altarelli-Parisi kernels for v ^ 0, 1 



PM = N C 
P q9 (v) = T R 



2 



4 + 2^(1 



v 2 + (1 - v 



(5.11) 



This result follows straightforwardly from the fact that in the limit considered 
we have y iS — > y^v and similar behaviour for other kinematical invariants. 
Using the phase-space integrals presented in appendix C we obtain the IR 
divergent contribution in the three-jet region of the four-parton processes, 



r^ = 2f°a s 2 JdPS(3)T bbg (y 13 ,y 23 ) 



2/12 



^y-c F +- e {c F -^ J {yi2+2rb)(3i2 



+N C 



log C12 



~3 Tr 




(N F -l)- + \ogr b 



(5.12) 



where, for consistency, we have also included a logr^ term coming from the 
integration of four massive quark transition probabilities. From eq. (4.7), 
eq. (5.12) we can easily see that, as expected, the IR divergences cancel be- 
tween three-parton and four-parton contributions rendering the final answer 
completely finite. In fact, one can see that this cancellation occurs separately 
for the different groups of diagrams defined in fig. 4. 
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6 Results and applications: Rf, D% e and m^mz) from data at the 
Z peak 



Since a big part of the calculation has been done numerically, it is important 
to have some checks of it. We have checked our transition probabilities for 
four-parton final states in the massless limit against the ones calculated by 
Ellis, Ross and Terrano (ERT) [69]. The massless limit cannot be taken di- 
rectly in the three-parton transition probabilities, since, as commented before, 
some collinear poles in e appear in our calculation as logarithms of the heavy 
quark mass. As seen before, we have checked that all the IR divergences cancel 
between three-parton and four-parton contributions in the massive case. To 
check the performance of the numerical procedure we integrated the massless 
amplitudes of ERT and obtained the known results for the functions In 
addition our four-parton amplitudes have been checked in the case of massive 
quarks, in four dimensions, by comparing their contribution to four-jet pro- 
cesses to the known results [33,76]. Finally, we have checked, independently 
for each class of diagrams with different colour factors, that the final result 
obtained with massive quarks reduces to the massless result in the limit of 
very small masses. 

The last check is the main check of our calculation. We have checked numer- 
ically that in the limit of M b -» we recover the values already known for 
the functions A^°\y c ) and A^\y c ) in the different algorithms considered. We 
did so by computing the functions Hy(y c ,rt)) and Ha^cH) for several small 
values of and then extrapolating the results for r& — > 0. This check is not 
trivial at all due to the already emphasized difference in the IR structure of 
the massless and massive cases [8]. 

Finally we have compared numerically some partial results with the results 
obtained in [38-41], where a completely different approach has been used to 
cancel IR divergences, and have found good agreement. 

In tables 9, 10, 11 and 12 we present the results for the functions Ay* and 
By^ A /Ay) for the different algorithms and different values of the pole mass of 

the 6-quark. For A$ we have presented the difference A$ — Ay\ which is 
entirely due to the triangle diagrams VI 2 in fig. 2. This difference is in general 
small and, as we discuss later, for the observables we are interested in, their 
contribution is suppressed. 

Using these functions or their derivatives we will study two interesting ob- 
servables that can be used to extract m\>{mz) from data at the Z peak. The 
observable proposed some time ago to measure the bottom-quark mass at the 
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Table 9 

The next-to-leading order functions (see eqs. (2.2, 2.4, 2.5)) for the E algorithm. 



Vc 


A 1 ) 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


165.014(156) 


874.1 ± 1.4 


467.5 ±0.7 


245.7 ±0.4 


0.02 


153.145(85) 


411.8 ± 1.0 


288.9 ±0.5 


193.8 ±0.3 


0.03 


128.583(64) 


221.4 ± 1.0 


192.4 ±0.5 


149.0 ± 0.3 


0.04 


107.624(53) 


136.1 ± 1.0 


130.2 ±0.5 


114.2 ±0.3 


0.05 


90.606(46) 


99.0 ± 1.1 


92.9 ±0.6 


87.6 ±0.3 


0.06 


76.839(40) 


78.2 ± 1.2 


72.5 ±0.6 


68.6 ±0.4 


0.07 


65.565(36) 


65.3 ± 1.2 


60.7 ±0.7 


56.4 ±0.4 


0.08 


56.198(32) 


57.4 ± 1.4 


52.8 ±0.8 


48.6 ± 0.4 


0.09 


48.344(30) 


50.4 ± 1.4 


47.0 ±0.7 


43.0 ±0.5 


0.10 


41.701(27) 


45.9 ± 1.8 


42.1 ±0.9 


38.9 ±0.5 


Vc 


A 1 ) a( l > 
A A V 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


0.97579(21) 


868.6 ± 1.4 


460.6 ±0.7 


239.7 ±0.4 


0.02 


0.89814(18) 


406.9 ± 1.0 


283.6 ±0.5 


188.2 ±0.3 


0.03 


0.82866(15) 


217.7 ±0.9 


188.2 ±0.5 


144.0 ± 0.3 


0.04 


0.76466(13) 


133.0 ± 1.0 


126.5 ±0.5 


109.8 ±0.3 


0.05 


0.70563(12) 


95.8 ± 1.1 


89.5 ±0.6 


83.7 ±0.3 


0.06 


0.65020(11) 


75.3 ± 1.2 


69.5 ±0.6 


65.0 ±0.4 


0.07 


0.59832(11) 


62.5 ± 1.2 


57.6 ±0.6 


53.2 ±0.4 


0.08 


0.54944(10) 


54.1 ± 1.4 


49.7 ±0.8 


45.4 ±0.4 


0.09 


0.50324(9) 


47.6 ± 1.4 


44.0 ±0.7 


39.9 ±0.5 


0.10 


0.45975(9) 


42.7 ± 1.8 


39.4 ±0.9 


36.1 ±0.5 



Z resonance was the ratio [7] 



where and T q are the three-jet and the total decay widths of the Z-boson 
into a quark pair of flavour q in a given jet-clustering algorithm. More precisely, 
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Table 10 

The next-to-leading order functions (see eqs. (2.2, 2.4, 2.5)) for the EM algorithm. 



Vc 


A 1 ) 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


139.008(152) 


-4.5 ± 1.0 


1.6 ±0.5 


4.4 ±0.3 


0.02 


131.434(82) 


-21.1 ±0.8 


-15.2 ±0.4 


-11.8 ±0.2 


0.03 


109.542(62) 


-19.4 ±0.8 


-16.0 ±0.4 


-13.3 ±0.2 


0.04 


90.576(51) 


-16.1 ±0.9 


-14.5 ±0.5 


-13.2 ±0.3 


0.05 


75.182(44) 


-13.7 ± 1.0 


-12.3 ±0.5 


-11.9 ±0.3 


0.06 


62.773(38) 


-11.6 ±1.1 


-11.0 ±0.6 


-10.7 ±0.3 


0.07 


52.685(35) 


-10.4 ± 1.2 


-10.0 ±0.6 


-10.0 ±0.4 


0.08 


44.354(31) 


-9.9 ± 1.3 


-9.4 ±0.7 


-9.2 ±0.4 


0.09 


37.443(28) 


-10.2 ± 1.5 


-9.0 ±0.7 


-8.9 ±0.4 


0.10 


31.650(26) 


-8.8 ± 1.6 


-8.2 ±0.9 


-8.3 ±0.5 


Vc 


A 1 ) a( l > 
A A V 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


0.97579(21) 


-4.8 ± 1.0 


0.7 ±0.5 


2.8 ±0.3 


0.02 


0.89814(18) 


-22.4 ±0.8 


-16.6 ±0.4 


-13.4 ±0.2 


0.03 


0.82866(15) 


-20.8 ±0.8 


-17.7 ±0.4 


-15.3 ±0.2 


0.04 


0.76466(13) 


-17.9 ±0.9 


-16.3 ±0.5 


-15.1 ±0.3 


0.05 


0.70563(12) 


-15.6 ± 1.0 


-14.4 ±0.5 


-13.8 ±0.3 


0.06 


0.65020(11) 


-13.6 ± 1.1 


-12.8 ±0.6 


-12.6 ±0.3 


0.07 


0.59832(11) 


-12.6 ± 1.2 


-11.9 ±0.6 


-11.6 ±0.4 


0.08 


0.54944(10) 


-12.1 ± 1.3 


-11.0 ±0.7 


-11.0 ±0.4 


0.09 


0.50324(9) 


-12.8 ± 1.5 


-10.6 ±0.7 


-10.4 ±0.4 


0.10 


0.45975(9) 


-11.8 ± 1.7 


-9.9 ±0.9 


-9.9 ±0.5 



the measured quantity is [14] 



3 ~ ri(y c )/r* 



where now we normalize on the sum over all light flavours (£ = u, d, s). To a 
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Table 11 

The next-to-leading order functions (see eqs. (2.2, 2.4, 2.5)) for the JADE algorithm. 



Vc 


A 1 ) 

A v 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


22.352(146) 


638.1 ±6.3 


533.2 ±3.2 


454.0 ± 1.8 


0.02 


71.356(78) 


27.7 ± 1.4 


28.1 ±0.7 


25.4 ±0.4 


0.03 


72.138(58) 


-4.7 ± 1.2 


-1.2 ±0.6 


-0.1 ±0.4 


0.04 


65.303(48) 


-12.7 ± 1.2 


-9.8 ±0.6 


-8.3 ±0.4 


0.05 


57.366(41) 


-15.1 ± 1.3 


-12.6 ±0.6 


-11.4 ±0.4 


0.06 


49.904(36) 


-16.7 ± 1.3 


-13.8 ±0.7 


-12.7 ±0.4 


0.07 


43.254(32) 


-16.7 ± 1.4 


-14.0 ±0.7 


-13.1 ±0.4 


0.08 


37.438(28) 


-16.5 ± 1.6 


-14.7 ±0.8 


-13.5 ±0.5 


0.09 


32.357(26) 


-16.3 ± 1.6 


-13.9 ±0.8 


-13.1 ±0.5 


0.10 


27.961(24) 


-16.4 ± 1.8 


-13.7 ±0.9 


-13.0 ±0.6 


Vc 


AW AW 

A A V 


B A ) / A V > 






M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


0.97579(21) 


661.7 ±6.3 


552.2 ±3.1 


468.1 ± 1.8 


0.02 


0.89814(18) 


29.8 ± 1.4 


29.2 ±0.7 


25.7 ±0.4 


0.03 


0.82866(15) 


-4.8 ± 1.2 


-1.4 ±0.6 


-0.6 ±0.4 


0.04 


0.76466(13) 


-13.2 ± 1.2 


-10.7 ±0.6 


-9.3 ±0.3 


0.05 


0.70563(12) 


-16.1 ± 1.3 


-13.6 ±0.6 


-12.5 ±0.4 


0.06 


0.65020(11) 


-18.0 ± 1.3 


-14.8 ±0.7 


-13.8 ±0.4 


0.07 


0.59832(11) 


-18.1 ± 1.4 


-15.3 ±0.7 


-14.3 ±0.4 


0.08 


0.54944(10) 


-18.6 ± 1.5 


-15.6 ±0.8 


-14.7 ±0.5 


0.09 


0.50324(9) 


-18.5 ± 1.6 


-15.7 ±0.8 


-14.4 ±0.5 


0.10 


0.45975(9) 


-19.4 ± 1.8 


-15.1 ±0.9 


-14.2 ±0.6 



good approximation both observables are related through 



r>M _ r>bd , a s A A Ay ( QA b _ Y^i=uds 9 A, \ / fi „\ 

3 3 - A (0) + [g Vb + g% E^uM, + 9%) ) ■ 1 J 

The values for A$ — Ay* can be taken from the tables. The extra term is 
mainly due to the fact that, contrary to R^ d , even in the M b — > limit the 
cancellation of triangle diagrams [68] is not complete in R^ e . The difference is 
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Table 12 

The next-to-leading order functions (see eqs. (2.2, 2.4, 2.5)) for the DURHAM 
algorithm. 



yc 










M h = 3GeV 


Mh = 4GeV 


Mh = 5GeV 


0.01 


38 229(56) 


14.7 ± 2.3 


12.6 ± 1.2 


q « _i_ n 7 


0.02 


34 531(28) 


-2.2 ± 1.7 


—2.5 ± 0.8 


-3.3 ± 0.5 


0.03 


28 488(20) 


-5.6 ± 1.6 


—4.5 ± 0.8 


—4.7 ± 0.5 


0.04 


23.521(16) 






t: . O _1_ VJ . O 


0.05 


19 574(13) 


-6.5 ± 1.8 


-4.9 ±0.9 


-4.7 ±0.6 


0.06 


16.425(11) 


-5.4 ±2.0 


-4.3 ±1.0 


-4.0 ±0.6 


0.07 


13 878(10) 


-3.9 ±2.1 


-3.6 ±1.1 


-4.0 ±0.7 


0.08 


11.786(9) 


-8.4 ± 2.4 


-4.8 ±1.2 


-4.1 ±0.7 


0.09 


10.045(8) 


-3.2 ±2.6 


-3.0 ±1.4 


-3.2 ±0.8 


0.10 


8.585(8) 


-6.1 ± 2.8 


-4.1 ± 1.5 


-3.3 ±0.9 


t/C 


A A 








M b = 3GeV 


M b = 4GeV 


M b = 5GeV 


0.01 


0.82567(22) 


18.5 ±2.3 


15.8 ±1.1 


12.1 ±0.6 


0.02 


70249(11) 


-0.8 ± 1.6 


-1.1 ±0.8 


-2.5 ±0.5 


0.03 


61124(9) 


-5.5 ± 1.6 


-3.7 ±0.8 


-4.2 ±0.5 


0.04 


53784(8) 


-8.1 ± 1.7 


-5.3 ±0.9 


-5.0 ±0.5 


0.05 


0.47625(7) 


-8.3 ± 1.8 


-5.5 ±0.9 


-5.4 ±0.6 


0.06 


0.42318(6) 


-7.7 ±2.0 


-4.9 ±1.0 


-4.4 ±0.6 


0.07 


0.37687(6) 


-6.8 ±2.1 


-4.5 ±1.1 


-4.5 ±0.6 


0.08 


0.33588(6) 


-8.4 ±2.4 


-5.9 ±1.2 


-3.9 ±0.7 


0.09 


0.29936(5) 


-5.5 ±2.6 


-4.7 ± 1.3 


-3.8 ±0.8 


0.10 


0.26654(5) 


-9.3 ±2.8 


-4.5 ± 1.5 


-3.7 ±0.9 



small and rather independent of the 6-quark mass. In the DURHAM algorithm 
and at y c = 0.02 it gives a contribution of ±0.002. In the other three algorithms 
considered this contribution is smaller than in DURHAM, since it is suppressed 
by a larger leading order function. Results for in the DURHAM 
algorithm have already been presented in [8,12]. Here we present B!§ for the 
four clustering algorithms discussed in this paper. 

This observable does not allow for a simultaneous analysis of a s and M b be- 
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cause the results for the different values of y c are correlated. To be able to 
make a fit for different values of y c we need a differential distribution. One can 
define the following ratio of differential two-jet rates 

D u = EM/c + Ay c /2) - I*(y c - Ay c /2)}/T b 
2 " [T%{y c + Ay c /2) - T%{y c - Ay c /2)]/r* " j 



For numerical results we fix Ay c = .005. The two-jet rate at 0(a 2 s ) is calculated 
from the three- and four-jet fractions through the identity eq. (2.1). As before, 
there is a small difference between D% e and the D\f ratio which is even smaller 
than the difference between R b £ and R\f. 

It is important to note that because the particular normalization we have used 
in the definition of these two observables most of the electroweak corrections 
cancel. Therefore, for our estimates it is accurate enough to consider tree-level 
values of gy and In addition, if they are expressed in terms of the effective 
weak mixing angle, most of the weak radiative corrections are also correctly 
taken into account [77]. 

From the definitions above, eq. (2.1), eq. (2.2), the values of the different 
functions given in the tables and using the known expression for T b (for mass 
effects at order a s see, for instance [7,56]) we obtain the values of the two 
observables as a function of y c and the quark mass. 

We use the following expansion in a s for the two observables considered 



Rf = 1 + ^ a (y c ) + r b L(r b , y c ) + ^ h(r b , y c )) , 
it y 7i J 

D bi = l + ?M c (y c ) + r b (d (r b , y c ) + ^ dfa, Vc)) ■ (6.5) 

7T y 7T j 

Note that the exact dependence on the heavy quark mass is kept in the func- 
tions bo, bi, d and d\, but for later convenience the leading dependence on 
r b has been factorized out. The a and c functions come exclusively from the 
triangle diagram and as commented before give a small contribution. 

Using the known relationship between the perturbative pole mass and the MS 
running mass [78], 



2a,(/i) f4 Mb 

TT \ 3 fl 2 



(6.6) 



we can re-express the same equations in terms of the running mass m b (/i). 
Then, keeping only terms of order 0(a s ) we obtain 
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a s (/x) _ / a s (/i) 
1 H an !/, + r h(u) On n, w r H 



^3 =H a (y c ) +r 6 (/i) 6 (r 6 ,y c )H h(r b , y c , fi) ] , 

7T y 7T y 

D* = l + ?M Co(yc) + fb(// ) yc ) + J l(f6> y a ^ , (6 .7) 

7T V 7T / 



where r b (/z) = m\{p)jm\ and 



4 _ /i 



2 



bi{r b ,y c ,fi) =bi(r b ,y c ) + 2b (r b ,y c ) [-- log r 6 + log- 2| . 

\ O 777- 5 



,2 



di(r b ,y c ,fi) = di(r b ,y c ) + 2d (r b ,y c ) - -logr 6 + log— ] . (6.8) 

\ O 777- r 



4 , ,„„ ^ 



The connection between pole and running masses is known up to order a 2 s1 
however, consistency of our pure perturbative calculation requires we use only 
the expression above. Although at the perturbative level both expressions, 
eq. (6.5) and eq. (6.7), are equivalent, they give different answers since different 
higher order contributions are neglected. The spread of the results gives an 
estimate of the size of higher order corrections. 

We have performed simple fits to the functions bo and do, describing the leading 
order behaviour of the and observables respectively. Also the and c$ 
contributions have been parametrized in terms of simple functions. The pairs of 
functions (bi,di) and (6i,Ji) give the NLO heavy quark mass corrections when 
a description in terms of the pole mass, eq. (6.5), or in terms of the running 
mass, eq. (6.7), is used. Although, by using the relationship in eq. (6.8). one 
can pass from one set of functions to the other, we have also fitted them 
independently. A Fortran code containing these fits can be obtained from the 
authors on request. 

In figs. 8 and 9 we present the results for the two observables studied, F$f and 
D\ , in the four jet-clustering algorithms considered. In all cases we plot the 
NLO results written either in terms of the pole mass, M b = 4.6 GeV [79,80], 
or in terms of the running quark mass at mz, m b {mz) = 2.83 GeV. The 
renormalization scale is fixed to /i — m z and a s {m z ) = 0.118. For comparison 
we also show R b ^ and Dtf at LO when the value of the pole mass, M b , or the 
running mass at mz, m b {mz), is used for the quark mass. 

Note the different behaviour of the different algorithms. In particular the E 
algorithm. As already discussed in [7], in this algorithm the shift in the res- 
olution parameter produced by the quark mass makes the mass corrections 
positive while from kinematical arguments one would expect a negative effect, 
since massive quarks radiate less gluons than massless quarks. Furthermore, 
the NLO corrections are very large in the E algorithm, see for instance fig. 10, 
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Fig. 8. The observable as a function of y c at the NLO for the four algorithms 
studied. The dotted lines give the observable computed at the NLO in terms of 
the pole mass M& = 4.6 GeV, while the dashed lines correspond to the use of the 
running mass mbimz) = 2.83 GeV instead. In both cases the renormalization scale 
is fixed at fx = mz, and a s (mz) = 0.118. For comparison we also plot the LO results 
for Mb = 4.6 GeV (lower solid lines) and mb(mz) = 2.83 GeV (upper solid lines). 

where we compare the following ratio 

size of NLO = NLO(m b ( mz )) -LQ(m b (m z) ) 

NLO(m b (m z )) - 1 ' 1 ; 



for i?3 in the four algorithms considered. The NLO correction in the E algo- 
rithm can be as large as 90 % of the LO result, a fact that probably indicates 
that it is difficult to give an accurate QCD prediction for it. For the JADE 
algorithm the NLO correction written in terms of the pole mass starts to be 
large for y c < 0.02. Note, however that the NLO correction written in terms of 
the running mass is still kept in a reasonable range in this region. DURHAM, 
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Fig. 9. Same as in fig. 8 but for the observable Dtf 
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Fig. 10. Comparative study of the relative size of NLO corrections in the four jet 
algorithms considered, see eq. (6.9). 
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in contrast, is the algorithm that presents a better behaviour for relatively low 
values of y c , while keeping NLO corrections in a reasonable range. 
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Fig. 11. The observable as a function of the renormalization scale [i for a fixed 
value of y c for the JADE and DURHAM algorithms. Labels are as in fig. 8. 

The theoretical prediction for the observables studied contains a residual de- 
pendence on the renormalization scale ji: when written in terms of the pole 
mass it only comes from the /z-dependence in a s (n), when written in terms 
of the running mass it comes from both a s (fi) and the incomplete cancella- 
tion of the /x-dependences between m^jj) and the logs of /i which appear in 
the MS expression. To give an idea of the uncertainties introduced by this we 
plot in fig. 11 the observable Rf as a function of fi for a fixed value of y c . 
Here we only present plots for the JADE and DURHAM algorithms, which, 
as commented above, have a better behaviour. We use the following one-loop 
evolution equations 

a(/i) = , mM = m b (m z ) , (6.10) 



where a(/x) = a s (/x)/7r, K — 1 + a(mz)A) log(/i 2 / m l) with f3 = (11 — 
2/3A^)/4, 70 = 1 and Np = 5 the number of active flavours, to connect 
the running parameters at different scales. 

For a given value of i?| we can solve eq. (6.5) (or eq. (6.7)) with respect to 
the quark mass. The result, shown in fig. 12 in the JADE and DURHAM 
algorithms for a fixed value of R^ e , depends on which equation was used and 
has a residual dependence on the renormalization scale jj,. The curves in fig. 12 
are obtained in the following way: first from eq. (6.7) we directly obtain for 
an arbitrary value of fi between mz and m^/10 a value for the bottom-quark 
running mass at that scale, mj,(/i), and then using eq. (6.10) we get a value 
for it at the Z-scale, m^mz)- Second, using eq. (6.5) we extract, also for 
an arbitrary value of \x between m z and m^/10, a value for the pole mass, 
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Fig. 12. Extracted value of mb(mz) from a fixed value of i?(f using either pole mass 
expressions (NLO-M^) or running mass expressions (NLO-mb(m^)) as explained in 
the text. Solid lines are obtained by using one- loop running to connect the result 
at different scales. Dashed lines use two-loop renormalization group equations. 

Mfe. Then we use eq. (6.6) at n — Mb and again eq. (6.10) to perform the 
evolution from \i = Mb to /i — mz and finally get a value for m^mz)- The 
two procedures give a different answer since different higher orders have been 
neglected in the intermediate steeps. The maximum spread of the two results 
can be interpreted as an estimate of the size of higher order corrections, that 
is, of the theoretical error in the extraction of mb(mz) from the experimental 
measurement of B!§. For the taken values of i?(f the spread one gets in rrib{mz) 
in both algorithms is roughly a little bit less than ±200 MeV. However, the 
spread of the results is strongly dependent on y c in the JADE algorithm while 
in DURHAM it is almost independent of it. A fact that, once more, shows the 
good theoretical behaviour of the DURHAM algorithm. 

Although, our observables are formally order 0(a s ) and therefore compati- 
ble with the use of one-loop renormalization group equations to connect the 
running parameters at different scales, as a consistency check of the result 
we have also repeated the analysis by using the two-loop running evolution 
equations [81] 



a(/i) = 



a(m z ) 



tf + a(m z )&i (L + a{m z )h l - L ) ' 
1 + a{m z ) gi 



m b {m z ) = m b {n) K 9 ' 

1 + a(/i) ci 

where L = log if and c x = g x - b x g^ with b x = fa/fa, gi = ji//3 and 



(6.11) 
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This corresponds to the dashed lines in Fig 12. In this case the spread of 
the results is enlarged to roughly ±280 MeV mainly due to the change in the 
NLO-M fe curve. If only the running mass expression, NLO-m^m^), is used the 
result is more stable when passing from 1-loop to 2-loop evolution equations. 



7 Conclusions 

In this paper we presented a detailed description of the next-to-leading order 
QCD calculation of the cross-section for heavy quark three-jet production in 
e + e~ annihilation. Complete quark mass effects have been taken into account. 
As a phenomenological application, we consider two different observables by 
using four different jet-clustering algorithms. The size of unknown higher order 
corrections is estimated by studying the residual renormalization scale depen- 
dence and the uncertainty that appears when the theoretical predictions are 
written in terms of the running mass or the pole mass of the produced heavy 
quark. In particular, we have carefully studied how to extract a value of the 
bottom-quark mass from the experimental measurement of these observables 
defined in several jet-clustering algorithms. These results have already been 
used in the measurement of the 6-quark mass at the Z peak as well as in the 
test of the flavour independence of QCD. 
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A One-loop integrals contributing to three-parton final states 

We consider here the basic one-loop integrals needed to compute the virtual 
corrections to the three-parton decay q — > P1+P2+P3 with all external particles 
on-shell, p\ = p\ = M£ and p\ = 0. By using the Passarino-Veltman [74] 
procedure we reduce any vector or tensor one-loop integral to a combination of 
different scalar one-loop integrals. Furthermore, the one-loop corrected three- 
parton amplitudes contribute at second order in the strong coupling constant 
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only via their interference with the tree-level amplitudes and, therefore, only 
the real parts of scalar one-loop integrals [71] are relevant for our calculation. 
Since one- and two-point functions are simple and well-known we skip their 
presentation. In this appendix we concentrate on the three- and four-point 
functions, especially on those with infrared divergences. Following results for 
the one- loop functions are understood to refer just to their real part. Details 
of the calculation have been presented in [11]. 



A . 1 Three-point functions 

We define the following set of three-propagator one loop integrals 



16tt 2 



s~im ( \ 4l—D f 

CUU/13)=/i J (2n)° k 2 [(k + p n ) 2 - M 2 ][(k - p 2 ) 2 - M 2 } ' 
f d D k 1 

C02(y 13 )=/i 4 D J (2?r) D pp + pis) 2 _ M 2] p + pi) 2 _ M 2] » 



% r d D k 1 

—CQZ{y^)=^ D J {2n) D k 2 {k _ pzm+pi) 2_ M 2 ] . (A- 1 ) 

d D k 1 



lQ n 2 C0A (y^) ^ ° J (2n) D [k 2 - M*][{k + qf - M 2 ][{k + p 12 Y - M ; 2 1 ' 
1 ; C05(y 12 )=^- D 1 ' !L ' L 1 



21 ' 



16tt 2 vy±z/ ^ J (2tt) d k 2 [{k + Pi y - M*][{k - P2 y - mi 

where Pij = Pi + Pj and = 2(pj • Pj)/s with s = q 2 . In general these 
functions depend on all the masses in the propagators, on the difference of 
external momenta squared and on the square of the external momenta itself. 
Since we have fixed masses and we have imposed the on-shell condition the 
only remaining relevant arguments of these functions are the two-momenta 
invariants y iS and y 12 , respectively. 

For a general scalar three-point function free of infrared singularities the loop 
integral can be performed in D = 4 dimensions and the result [71] can be 
expressed as a sum over twelve dilogarithms (Spence functions). This is the 
case for the three-point integrals C01, C02 and C04. However, after some 
algebra, the real parts of the C02 and C04 functions take a simpler form 



C02(y 13 ) ' 



syi3 



1 , 2 t • / \ 71-2 

-log y + Li 2 (y ) - — 
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C04(y 12 ) = — 
where 

_ 1-/? 
1 -/3i2 

Cl2 



(2/12 + 2r b ) 2s 



log 2 C - log 2 C12 



/3 = v / r^4^, 



Pl2 = \l- 



Ar h 



I + A2' V* yi2 + 2r fe ' 

and |/o = r b / (r& + y 13 ). For the C01 function we get 
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x 2 = 



Z2 



a 2/13 
1 - ay/\ ' 

2/13 



^3 = 



yi3 



2/13 + n 

with A = (1 - 2/13) 2 - Ar b and a = (1 - y 13 + v / A)/2. 



(A.3) 



(A.4) 



For the C03 function, which has a double infrared pole, the result we obtain 
for the real part, after integration in D = A — 2e dimensions, is 



C03(y 13 ) = s- l - e fi 



l-e..2e (47T) 



T(l - e) 2/13 



X 



1 1 1 Tb , 1 1 2 r b 1,2 r • / \ 

2^ + 21 bg ^ + 4 bg S " 2 bg y ° " L ' 2(2/0) 



(A.5) 
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Notice that for convenience we have factorized a (47r) e /T(l — e) coefficient 
because this constant appears also in the four-body phase-space, see eq. (B.5). 

The C05 function has a simple infrared pole and already appeared in the NLO 
calculation of two-jet production, although in a different kinematical regime. 
Here we have 
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p r(l-e)(y 12 + 2r 6 )/3 12 
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Q - logr^ logc 12 - 2L(y 12 ) 



(A.6) 



where 

7T 2 1 2 - 

L(yi 2 ) = Li 2 (c 12 ) + y + log(l - ci 2 ) logci 2 - - log 2 c i2 . (A. 7) 

A . 2 Four-point functions 



Only two types of box integrals appear in the virtual corrections to the three- 
parton decay. We define 



D05(y 13 ,y 12 ) = (A.i 
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d D k 
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and 



D06(y 13 ,y 23 ) = (A.9) 



16tt 2 

,,4-d f d k 
* J (2n) D 



(2n) D k\k + p 3 ) 2 [(A; + P13) 2 - M 2 ] [(/c - p 2 ) 2 - M 2 } ' 
where p 13 = p 1 + p 3 and = 2(pj • Pj)/s with s = q 2 . 

Both four-point functions are IR divergent. The integral D05 has a simple 
infrared pole, while D06, that involves a three-gluon vertex, exhibits a double 
infrared pole. Furthermore, their divergent behaviour is closely related to the 
infrared structure of the C03 and C05 functions defined before in a very simple 
way [82,83] 



s ■ D05(y 13 , y 12 ) = —C05(y 12 ) + finite terms , 

J/13 

s ■ D06(y 13 , y 23 ) = —C03(y 23 ) + —C03(y 13 ) + finite terms . (A. 10) 

1/13 1/23 



The box integral, D05, is also met in one-loop electroweak calculations and it 
has been calculated [84] by using a photon mass infrared regulator. We just 
quote the result in dimensional regularization 
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Finally, we have the following result for the real part of the D06 function [11] 



£06(y 13 , y 23 ) = *- 2 -y e wr^ 1 (A - 12) 

r(l - e) ywy 23 
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B Phase-space in D = 4 — 2e dimensions 



The phase-space for n-particles in the final state in arbitrary space-time di- 
mensions, D, has the following general form [72,73] 



dPS(n) = (2*) n CTi" 2 ^ D [9 - E ft ] (B-l) 

i=l,n V / * V i=l,n 

1 — 1,71 V / \ 1 — 1,71 



Let's consider the decay into three particles, q P\+p 2 +p 3 , where particles 
1 and 2 share the same mass, p\= p\ = M 2 , and particle 3 is massless, p\ = 0. 
In terms of the two-momenta invariants y i3 = 2(p 1 -p 3 ) / s and y 23 = 2{p 2 -p 3 ) / s, 
with s = q 2 , we get (D = 4 — 2e) 



where the function h p , which defines the phase-space boundary, has the form 
hp = 1/131/23(1 - 1/13 - I/23) - r b (y 13 + y 23 f . (B.3) 

For the case of the decay into two massive and two massless particles, q 
Pi + P2 + Ps + P4 with p\ = p\ = M 2 and p\ = p\ = 0, and when the two 
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massless particles can become collinear, it is convenient to write the four-body 
phase-space as a quasi three-body decay 



q -> P34 + Pi + P2 
^ P3 + Pa ■ 

In the cm. frame of particles 3 and 4 the four-momenta can be written as 



Pl = (£?!,..., 0,pi) , 

P2 = (-E2,---,P2sin/3,p 2 cos/3) , 

Ps = £3(1, . . . , sin 6 1 cos cos#) , 

p 4 = £4(1, . . . , — sin 6 cos 9' , — cos #) , 

where the dots in p 3 and p 4 indicate .D — 3 unspecified, equal and opposite 
angles (in D dimensions) and D — 3 zeros in pi and pi- We will refer to this 
as the "3-4 system" [69]. 



In terms of the following invariants 
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where p 34 = p 3 +p 4 , energies and three-momenta become in the "3-4 system" 
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Defining t> = (1 — cos#)/2, the D-dimensional phase-space in this system is 
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where a 1/2 statistical factor has been included, Ngi is a normalization factor 

r(l-2e) 
r 2 (l-e) ' 



Ny = J dff sin" 2 * ff = 2 2 -7T ^ i _ 2 e £) , (B.6) 
o 



and the function 

h 3 4 = 2/1342/234(1 - 2/134 - 2/234) - n(y 134 + 2/234) 2 

- (l - Ar b - 2(1 - 2r b )(y lu + y 2M ) + 2/ 2 34 + 3y lu y 23 4 + 2/234) 2/34 

+ 2(1 - 2r b - y lu - 2/234)2/34 - Vu » ( B - 7 ) 

defines the limits of the phase-space. 

For the integration of the parts of the four-parton transition probability con- 
taining only quark-gluon soft singularities it is convenient to use another 
parameterization of the phase-space. For example, if the invariant 2/13 — 
2(pip 3 )/s approaches zero, we choose the so-called "1-3 system" [69]. Intro- 
ducing variables 

2(pi-p 3 ) 2(p 2 -p 13 ) 2(p 4 -p 13 ) 

2/13 = , 2/213 = , 2/413 = , (B.8) 



with pi 3 = pi +p 3 . Energies and three- momenta in this system read 
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In this system, the four-body phase-space has the form 
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where now the angle 9 is the one defined by the three- momenta pi and P2, and 
the function that defines the limits of the phase-space is 



hi3 = (1/2131/413 + 4rf)(l - y 2 i 3 - 2/413) - 4rf 

-r b (l- 2(y 2 13 + 1/413) + 2/213 + %213l/413 + %413) 

- (l - Qr b + 8r^ - 2(1 - 3r b )(y 2 i 3 + 2/413) + 2/213 + 3 l/2i3l/4i3 + I/413) S/13 
+ (2 - 5r 6 - 2y 213 - 2y 4 i 3 ) y? 3 ~ 3/?3 • (B.10) 



C Infrared divergent phase-space integrals 

We consider now the tree-level four-parton decay q — > Pi + P2 + P3 + Pa- 
When at least one gluon is radiated from an external quark (or gluon) line, 
the propagator-factors, l/(pi -Pj), can generate soft and collinear divergences 
at the border of the four-body phase-space. In this appendix we present the 
basic results that appear when these factors are integrated in a thin slice at the 
border of the phase-space and we show how these basic phase-space integrals 
are related to the infrared divergent scalar one-loop integrals discussed in 
appendix A. 

C. 1 Integrals containing soft gluon divergences 

Phase-space of n + 1 particles can be written as the product of the n-body 
phase-space times the integral over the energy and the solid angle of the extra 
particle. In arbitrary D = 4 — 2e dimensions we have 



Suppose E 3 is the energy of a soft gluon. Assuming that E 3 < w^fs, where w 
is very small and defines an upper cut on the soft gluon energy. Then we have 
the following useful results 
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Notice that this last integral has the same divergent structure as the scalar one- 
loop C05 function defined in eq. (A. 6). The finite contribution, the F(y 1 4, 2/24) 
function, is rather involved. We write it in terms of the variables y± 3 and y 23 
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where h p is the function that defines the limits of the three-body phase-space, 
see eq. (B.3), and 



fi, 2 = (1 ± VT^)/a , 
z\ = exp 

z 2 = exp 



(C.7) 
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C.2 Integrals containing gluon-gluon collinear divergences 



For transition probabilities containing gluon-gluon collinear divergences we use 
the four-body phase-space representation of eq. (B.5). In the limit j/ 34 — > 
the function that defines the boundary of the four-body phase-space, /i 34 in 
eq. (B.7), reduces to the three-body phase-space function h p in eq. (B.3) and 
P34 — P3 + £>4 behaves as the momentum of a pseudo on-shell massless particle 
since — > 0. Therefore, in this limit, an effective three-body phase-space 
can be factorized in the same way we made for treating the soft singularities. 
After integrating over the two angular variables and over the 1/34 invariant for 
I/34 < w we get the following result 
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Notice that we get the same infrared poles as in the one-loop three-point 
function C03, eq. (A. 6), if we identify j/134 with y± 3 . 
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